/*=============================================================================
*
*   Copyright (C) 2020 All rights reserved.
*
*   Filename:		triangle_shape_function.hpp
*
*   Author: Wang Zhecheng - wangzhecheng@yeah.net
*
*   Date: 2020-03-22 22:57
*
*   Last Editors: Wang Zhecheng - wangzhecheng@yeah.net
*
*   Last modified:	2020-08-04 01:39
*
*   Description:
*
=============================================================================*/
#ifndef TRIANGLE_SHAPE_FUNCTION_H
#define TRIANGLE_SHAPE_FUNCTION_H

#include "../../types.hpp"

namespace geoxfem {
namespace geometry {

template <size_t dim, size_t degree>
struct triangular_shape_function {
  static constexpr _real r_5 = .5;
  template <typename Shape>
  static VectorXR getStdNByLocal(const VectorXR&, const Shape&) {
    throw "not implemented yet";
  }
  template <typename Shape>
  static MatrixXR getStddNdlxByLocal(const VectorXR&, const Shape&) {
    throw "not implemented yet";
  }
  template <typename Shape>
  static VectorXR getStdNByGlobal(const VectorXR&, const Shape&) {
    throw "not implemented yet";
  }
  template <typename Shape>
  static VectorXR global2local(const VectorXR&, const Shape&) {
    throw "not implemented yet";
  }
};
/*triangle| dimension: 2D, order: 1{{{*/

template <>
template <typename Shape>
VectorXR triangular_shape_function<2, 1>::getStdNByLocal(const VectorXR& lcoord,
                                                         const Shape&) {
  Vector3R    N;
  const _real xi  = lcoord(0);
  const _real eta = lcoord(1);
  N << xi, 1 - xi - eta, eta;
  return N;
}
template <>
template <typename Shape>
MatrixXR triangular_shape_function<2, 1>::getStddNdlxByLocal(const VectorXR&,
                                                             const Shape&) {
  MatrixXR dNdlx(2, 3);
  dNdlx << 1, -1, 0,  //
    0, -1, 1;         //
  return dNdlx;
}
template <>
template <typename Shape>
VectorXR
triangular_shape_function<2, 1>::getStdNByGlobal(const VectorXR& gcoord,
                                                 const Shape&    shp) {
  Matrix3R    m;
  Vector3R    v;
  const _real x1 = shp.at(0).x();
  const _real y1 = shp.at(0).y();
  const _real x2 = shp.at(1).x();
  const _real y2 = shp.at(1).y();
  const _real x3 = shp.at(2).x();
  const _real y3 = shp.at(2).y();
  m << 1, 1, 1,  //
    x1, x2, x3,  //
    y1, y2, y3;
  v << 1, gcoord[0], gcoord[1];
  return m.inverse() * v;
}
template <>
template <typename Shape>
VectorXR triangular_shape_function<2, 1>::global2local(const VectorXR& gcoord,
                                                       const Shape&    shp) {
  Vector2R lcoord;
  MatrixXR refCoord(2, 3);
  refCoord << 1, 0, 0,  //
    0, 0, 1;
  VectorXR N = getStdNByGlobal(gcoord, shp);
  lcoord << refCoord * N;
  return lcoord;
}

/*}}}*/
/*triangle| dimension: 2D, order: 2{{{*/

template <>
template <typename Shape>
VectorXR triangular_shape_function<2, 2>::getStdNByLocal(const VectorXR& lcoord,
                                                         const Shape&) {
  VectorXR    N(6);
  const _real c1 = lcoord(0);
  const _real c2 = lcoord(1);
  const _real c3 = 1 - c1 - c2;

  N << (2. * c1 - 1.) * c1,  //
    4. * c3 * c1,            //
    (2. * c3 - 1.) * c3,     //
    4. * c2 * c3,            //
    (2. * c2 - 1.) * c2,     //
    4. * c1 * c2;

  return N;
}
template <>
template <typename Shape>
MatrixXR
triangular_shape_function<2, 2>::getStddNdlxByLocal(const VectorXR& lcoord,
                                                    const Shape&) {
  MatrixXR    dNdlx(2, 6);
  const _real c1 = lcoord(0);
  const _real c2 = lcoord(1);
  const _real c3 = 1 - c1 - c2;

  dNdlx << 4. * c1 - 1.,  //
    4. * (c3 - c1),       //
    -(4. * c3 - 1.),      //
    -4. * c2,             //
    0.,                   //
    4. * c2,              //
    0.,                   //
    -4. * c1,             //
    -(4. * c3 - 1.),      //
    4. * (c3 - c2),       //
    4. * c2 - 1.,         //
    4. * c1;
  return dNdlx;
}
template <>
template <typename Shape>
VectorXR
triangular_shape_function<2, 2>::getStdNByGlobal(const VectorXR& gcoord,
                                                 const Shape&    shp) {
  MatrixXR    m(6, 6);
  VectorXR    v(6);
  const _real x1 = shp.at(0).x();
  const _real y1 = shp.at(0).y();
  const _real x2 = shp.at(1).x();
  const _real y2 = shp.at(1).y();
  const _real x3 = shp.at(2).x();
  const _real y3 = shp.at(2).y();
  const _real x4 = shp.at(3).x();
  const _real y4 = shp.at(3).y();
  const _real x5 = shp.at(4).x();
  const _real y5 = shp.at(4).y();
  const _real x6 = shp.at(5).x();
  const _real y6 = shp.at(5).y();
  const _real x  = gcoord(0);
  const _real y  = gcoord(1);
  m << 1, 1, 1, 1, 1, 1,                                   //
    x1, x2, x3, x4, x5, x6,                                //
    y1, y2, y3, y4, y5, y6,                                //
    x1 * x1, x2 * x2, x3 * x3, x4 * x4, x5 * x5, x6 * x6,  //
    x1 * y1, x2 * y2, x3 * y3, x4 * y4, x5 * y5, x6 * y6,  //
    y1 * y1, y2 * y2, y3 * y3, y4 * y4, y5 * y5, y6 * y6;
  v << 1, x, y, x * x, x * y, y * y;
  return m.inverse() * v;
}
template <>
template <typename Shape>
VectorXR triangular_shape_function<2, 2>::global2local(const VectorXR& gcoord,
                                                       const Shape&    shp) {
  Vector2R lcoord;
  MatrixXR refCoord(2, 6);
  refCoord << 1, 0.5, 0, 0, 0, 0.5,  //
    0, 0, 0, 0.5, 1, 0.5;
  VectorXR N = getStdNByGlobal(gcoord, shp);
  lcoord << refCoord * N;
  return lcoord;
}

/*}}}*/
////////////////////////////////////////////////////////
////////////////////////////////////////////////////////

}  // namespace geometry
}  // namespace geoxfem

#endif /* TRIANGLE_SHAPE_FUNCTION_H */
